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Abstract. The gas-phase chemistry of dark clouds has been studied with a treatment of uncertainties caused both 
by errors in individual rate coefficients and uncertainties in physical conditions. Moreover, a sensitivity analysis 
has been employed to attempt to determine which reactions are most important in the chemistry of individual 
species. The degree of overlap between calculated errors in abundances and estimated observational errors has 
been used as an initial criterion for the goodness of the model and the determination of a best 'chemical' age of the 
source. For the well-studied sources L134N and TMC-1CP, best agreement is achieved at so-called "early times" 
of ~ 10 5 yr , in agreement with previous calculations but here put on a firmer statistical foundation. A more 
detailed criterion for agreement, which takes into account the degree of disagreement, is also proposed. Poorly 
understood but critical classes of reactions are delineated, especially reactions between ions and polar neutrals. 
Such reactions will have to be understood better before the chemistry can be made more secure. Nevertheless, 
the level of agreement is low enough to indicate that a static picture of physical conditions without consideration 
of interactions with grain surfaces is inappropriate for a complete understanding of the chemistry. 
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1. Introduction 

The chemistry of dark clouds, now known as qui- 
escent cores, has been studied for ove r thirty years 
(|Herbst k Klempererlll973 IWatsonlll973|) . Indeed, a sig- 
nificant fraction of the more than 130 gas-phase species 
detected via high-resolution spectroscopy in the inter- 
stellar and circumstellar media has been seen in such 
sources, whi ch include the well - studied clouds TMC1-CP 
and L134N llSmith et alJl2004t iTerzieva k Herbstlll99Sli . 
In an attempt to reproduce the abundances of these 
species, a large number of gas-phase chemical models have 
been developed that compute the variation of the gas- 
phase concentrations as a function of time. For many 
years, most of these models w ere based on the pseudo- 
time-depend e nt ap p roximation JpTasad k Huntresslll980t 
Leung et alJ Il984t lLanger et alJ Il984 iMillar k Neiadl 



19851) . in which the chemistry evolves under fixed and 



homogeneous conditions from some initial abundances, 
consisting of totally atomic material except for a high 
abundance of H2. In the last decade, more complex 
models, including surface chemistry l|Hasegawa et all 
Il992t iRuffle k Herbsd l2000|) . heterogeneity of the source 
l|BhirtaTu^t^t^iTT2^01 ) . and assorted effects of dynam- 



ics and turbulence l|Markwick et alj EoOO) , have been re- 
ported. 

In all of these models, the connections among gas- 
phase species are described by up to thousands of chem- 
ical reactions, with rates quantified by rate coefficients, 
most of which are poorly determined by experimental or 
theoretical means. Nevertheless, even the older pseudo- 
time-dependent models were at least partially success- 
ful in reproducing the exotic nature of the chemistry 
(molecular ions, radicals, and metastable isomers), the 
unsaturated nature of the more complex products (e.g., 
cyanopolyynes) , and the strong deuterium fractionation. 
They have been less successful, however, in detailed com- 
parisons with observations of fractional abundances and 
column densities for the large number of species observed 
in individual sources. In the last published comparison 
between a pseudo-time-dependent the ory and observatio n 
for the well-studied source TMC1-CP l|Smith et alJl2004|h 
it was found that at the time of best agreement (1 x 10 5 
yr), only about 2/3 of the observed species have calculated 
fractional abundances within an order-of-magnitude of the 
observed values. This result, obtained with the relatively 
new osu.2003 network and so-called "low- metal" elemental 
abundances, is actually worse than obtained with a pre- 
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abunda nces were reproduced to within one order of mag- 
nitude ijTerzieva fc Herbsllll998j) . The order-of-magnitude 
criterion used in both studies is a subjective one, however, 
and no comparison involving TMC1-CP has ever been 
done by taking into account the rate coefficient uncer- 
tainties in the computation of the theoretical abundances. 

Without any quantification of the random model er- 
ror, it is difficult to conclude if observed abundances can 
be reproduced or not by the model. Indeed, including the 
estimated uncertainties in rate coefficients is the only way 
to define those species with abundances not reproduced 
by chemical networks. With a rigorous quantitative ap- 
proach to uncertainty, moreover, one can begin to compre- 
hend the deficiencies of gas-phase pseudo-time-dependent 
calculations, and determine where improvement is most 
necessary. If, for example, hydrogen-rich species, which 
may be synthesized on grain surfaces and desorbed into 
the gas via non-thermal means, are the only species with 
abundances to be poorly reproduced, the most important 
correction would involve the inclusion of surface processes. 
If, on the other hand, smaller species such as O2 and H2O, 
with simple chemistry, cannot be reproduced within the 
estimated errors, then it is at least conceivable that a far 
more complex physical scenario is needed in which the for- 
mation and dynamics of the quiescent core play roles in 
the chemistry (Bergin & Melnick 2005). Whether or not 
to give up on the pseudo-time-dependent approximation 
then depends critically on our estimates of uncertainties, 
because there is no point in violating Occam's Razor if a 
more simple theory is as adequate as a more complex one. 

This p aper is the second ap plication of a method de- 
veloped in IWakelam et al.l ( 20051) to quantify the ran dom 
errors of chemical models. In IWakelam et ahl l|2005l) . we 
presented the method in great detail and showed some 
consequences for hot core chemistry (T> 100 K). Hot 
cores occur during an evolutionary stage of star forma- 
tion, a nd can have complex s tructures not well under- 
stood l|Nomura fc MilladboO^) . Dark clouds, which have 
not yet started to collapse in an appreciable manner, may, 
on the contrary, have relatively simple temperature and 
density structures. These objects are then better labora- 
tories to test the gas-phase chemical models usually used. 
In this paper, we apply the method to estimate the statis- 
tical errors as a function of time for dark cloud chemistry, 
which is occurring at an H2 density of ~ 10 4 cm -3 and a 
temperature < 20 K. A previous analysis of uncertainties 
for mainl y steady-state conditio ns was performed for dark 
clouds bv lVasvunin et all 1 20041) . although their error was 
defined in a different manner. 

The remainder of our paper is organized as follows. In 
Section 2, we briefly describe the chemical network and ap- 
proach to uncertainties used for this study. Some general 
results are presented in Section 3, where we also (i) con- 
sider differences between the two major networks, and (ii) 
make a com parison between our un certainty calculations 
and those of IVasvunin et af] l|2004t) . Section 4 contains a 
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Fig. 1. Mean abundance (< X >) (left panels) and error 
AlogX (right panels) for N and SO plotted as a function 
of the number of runs. Both quantities are normalized 
by the value obtained for the maximum number of runs 
(2500). The results are given for Model 2 at two times: 
10 4 yr on the lower panels and 10 6 yr on the upper panels. 



Table 1. Physical parameters for the models. 



Physical conditions 
10 K, 10 4 cm- 3 " 



Model 1 

Model 2 5-15 K, (0.5 - 1.5) x 10 4 cm" 3 " 

Model 3 5-15 K, (0.5 - 1.5) x 10 4 cm" 3 ", C/0=1.2 

a H 2 density. 



Table 2. Initial elemental abundances with respect to H2. 



Initial abundances 



He 


2.8 x 10" 1 


Fe+ 


6.0 x 10" 


a 


N 


4.28 x 10 -5 


Na+ 


4.0 x 10" 


9 


O 


3.52 x 10~ 4 


Mg+ 


1.4 x 10" 


8 


c+ 


1.46 x 10" 4 


P+ 


6.0 x 10" 


9 


S+ 


1.6 x 10" 7 


C1+ 


8.0 x 10" 


9 


Si+ 


1.6 x 10 -8 









Section 5 discusses the specific cases of O2 and H2O, and 
the cloud ages. We conclude the paper in Section 7. 

2. Chemical network and uncertainty method 

We used a time-independent physical model with 
the gas-pha s e che mical network osu.2003 1 reported by 
ISmith et alJ l|2004l) . This network contains standard gas- 
phase reactions (ion-neutral, neutral-neutral, dissociative 
recombination etc), with the addition of a significant num- 
ber of rapid radical- neutral processes. Except for H2 pro- 
duction, the grains are only important as sites of negative 
charge that recombine with positive ions. We considered 
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the uncertainties in all the classes of reactions except for 
ion recombination with negative grains. The uncertainties 
in the rate coefficients have recently been added to the 
osu.2003 network, based on t hose listed in the RATE99 
network llLe Teuff et al.ll2000R For the uncertainties not 
estimated at 10 K (according to the temperature range 
given in RATE99), we increased the uncertainty to a fac- 
tor of 2. 

The Monte Carlo method used here to include the 
rate-coefficient uncer tainties is described in detail by 
I Wakelam et alJ l)2005|) . Briefly, it consists of generating N 
new sets of rate coefficients by replacing each coefficient fcj 
by a random value consistent with its uncertainty factor 
Fi. We assume a normal distribution of log fc* with a stan- 
dard deviation Oi = log Fi. We run the model for each 
set j, which produces, for each species, N values of the 
abundance Xj(t) at a time t. The mean value of logX(i) 
gives us the "recommended" value while the dispersion 
of log Xj(t) around < logA(i) > determines the error 
due to kinetic data uncertainties. The error in the abun- 
dance A log A = ^(logA max — logA min ) is defined at a 
time t as the smallest interval [logA min , logX max ] that 
contains 95.4% of log A.,- values, which is equivalent to 
2r7 for symmetric Gaussian distributions. For example, if 
A log A = 1.0 and the calculated distribution is symmet- 
ric, then the 2a values of A max and A m ; n are one order of 
magnitude greater and less than X, respectively. Similarly, 
if AlogX = 0.5, the 2<r values are greater and less than 
X by a factor of 3.3, and if A log X = 0.3, they are greater 
and less than X by a factor of 2.0. 

We initially ran two models with different physical con- 
ditions. In the first model (Model 1), the temperature and 
the H2 density are held constant at values of 10 K and 
10 4 cm -3 , respectively; only the rate coefficients vary. For 
the second model (Model 2), we randomly chose the tem- 
perature and the density within a possible range of values. 
This uncertainty was adopted for two reasons: (i) the cloud 
temperature and density are usually derived from obser- 
vations (using approximate line excitation models) and an 
error in these values can usually be determined, and (ii) 
the modeled clouds are more likely to be inhomogeneous 
sources with contributions over the chosen ranges of tem- 
perature and density. We investigated the sensitivity of the 
results if an uncertainty of ±50% is considered for T and 
n(H2) around the typical values given for Model 1 (see 
Table ^ . Because temperature and density are physical 
parameters not characterized by a Gaussian distribution, 
we used a flat distribution instead for the rate coefficients. 
For both models, we used a cosmic-ray ionization rate of 
1.3 x 10~ 17 s _1 , a visual extinction of 10 and the low- 
metal elemental abundances, which are listed in Table [21 
A third model, with carbon-rich elemental abundances, is 
discussed later in the text. 

For Model 1 and Model 2, we performed 2000 and 2500 
different runs, respectively. In order to demonstrate con- 
vergence, we plot as examples in Fig. ^ the mean abun- 
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Fig. 2. Error A log A as a function of time for some 
species with Model 1. 
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Fig. 3. Error, A log A, as a function of the number of 
atoms per molecule at 10 6 yr for Model 1. 



dance < A > and its error Alog(A) for the species N 
and SO as functions of the number of runs at two dif- 
ferent times with Model 2. The plotted parameters are 
normalized by the values obtained for 2500 runs. For all 
the species, we noticed that the number of runs chosen is 
more than adequate. 

3. Calculated uncertainties in abundances 

3.1. General considerations 

Figure |2 presents the error A log A for some species as a 
function of time and for the physical conditions of Model 
1 (see Table . As already noticed for hot core chem- 
istry, the errors in the abundances increase with time for 
most of the molecules except the species that contain a 
dominant portion of an element, such as CO for carbon 
and N2 for nitrogen. For example, the errors for HCN and 
HC3N are 0.2 and 0.5 at 10 5 yr and increase to 0.3 and 
1.1, respectively, at 10 7 yr. The error also increases with 
th e complexity of the mol ecule, an effect previously noted 
bv IVasvunin et aTl l)2004[) . This effect is shown in Fig. [3] 
for the physical parameters of Model 1 at a time of 10 6 
yr. For molecules with 10 atoms, the figure shows that an 
error of 1.0, corresponding to a factor of 10, is in the mid- 
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Fig. 4. Density of probability of the NHJ abundance. The 
right plots represent the histograms of the abundance dis- 
tributions at 10 7 yr. The top row is for Model 1 and the 
bottom for Model 2. 
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Fig. 5. Error, AlogX, as a function of the number of 
atoms per molecule at 10 6 yr for Model 2. 



3.2. Variation of the temperature and density 

The variations of the gas temperature and density in 
Model 2 cause two major effects. First, they modify the 
distribution of the abundances at a given time from a 
Gaussian shape to an asymmetrical one, as can be seen by 
the example of NH4" in Fig. 21 The effect is not surprising 
since the rate equations are not linear with respect to den- 
sity or temperature. Indeed, the dependence on tempera- 
ture can be exponential for processes with a small barrier. 
The second effect is an increase in the error, which can 
be significant for some of the species, as shown in Fig. [5] 
which is to be compared with Fig. [3] The nitrogen-bearing 
species are particularly sensitive to the variation of tern- 



Table 3. Examples of errors computed using Model 1 
(AlogXi) and Model 2 (AlogX 2 ) at 10 5 yr. 



Species 


AlogXi 


A log X 2 


N+ 


0.36 


1.0 


HNO 


0.31 


0.91 


NH 3 


0.27 


0.88 


NH 2 CN 


0.37 


0.87 


SO 


0.64 


0.79 


H 2 


0.40 


0.42 


HC 3 N 


0.45 


0.56 


C 9 N 


0.43 


0.53 



that the variation in temperature is the cause because the 
nitrogen chemistry starts with the reaction N + + H2 — > 
NH + + H, which is assumed to be endothermic by 85 K in 
both the osu.2003 and RATE99 networks. The situation is 
more complex than can be contained in networks because 
it is likely that an estimated non-thermal fraction of H 2 
in its J=l ortho state of 10~ 3 drives the reaction and es- 
sentially converts all the atomic nitrogen ion into NH + 
(|Le Bourlodll99l . Nevertheless, with the adopted rate 
coefficient, the reaction to form NH + is so inefficient at 
temperatures under 10 K that it loses importance. Thus, 
temperatures lower than 10 K lead to much lower abun- 
dances for some N-bearing species. The result is a skew- 
ing of the abundance distributions to lower values, as can 
be seen for the case of NH|" in Fig. 0] The effect of the 
temperature and density variations is less important for 
complex molecules because the variations caused by the 
physical parameters are diluted by the large uncertainties 
due to the uncertainties in rate coefficients. For example, 
the uncertainties in the cyanopolyne abundances do not 
show a strong dependence on the physical parameters: the 
increase of A log X with Model 2 is less than a factor of 
2 for HCN and this number decreases with the increasing 
complexity of the cyanopolyne. 

3.3. Comparison between osu.2003 and RATE99 
databases 

The statistical uncertainties in the rate coefficients are not 
the only sources of error for gas-phase models. Most of the 
reactions have not been studied in the laboratory or via 
detailed theoretical considerations, so that approximate 
values for the rate coefficients must be used. Even studied 
reactions, if studied by more than one group, show large 
discrepancies in rates. Some of the problems in choosing 
proper rate coefficients can be seen in a comparison of 
the osu.2003 and RATE99 databases, which contain many 
reactions with different choices of rate coefficients. These 
differences arise from several specific sources: 

a) Different estimates for poorly understood reactions 
such as radical-neutral reactions and both neutral- 
neutral and ion- molecule radiative associations, 
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Fig. 6. CO, H2O and CS abundances as a function of 
time computed using the osu.2003 (dashed lines) and the 
RATE99 (solid lines) databases. The three curves for each 
molecule and database refer to the average value and the 
2a errors. 



Database (http://kinetics.nist.gov/index.php) or from 
competing measurements, 
c) Different approximations regarding the temperature 
dependence of ion-polar neutral reactions. 

The p roblem with radical- neutral reactions has been stud- 
ied bv I Smith" et all l|2004|) : the osu.2003 network includes 
estimates of rapid rates for a variety of such reactions that 
are greater than used in RATE99 as well as the previous 
network from the Ohio group. The last discrepancy derives 
from the fact that although experiments show that virtu- 
ally all of the small number of ion-polar neutral rea ctions 
studied possess an inverse tem perature dependence ([Rowel 
Il992t iRowe fc Rebrionlll992h this dependence is not eas- 
ily expressible in terms of the simple rat e expression used 
in the networks. Based on the work of iHerbst k. Leunel 
l)l986|) . the osu.2003 network uses an approximation de- 
rived from the locked dipole approximation for linear neu- 
trals and classical scaling ap proach for non-linear neu- 
trals l)Su &: Chesnavichlll982[) . Both of these approxima- 
tions lead to a temperature dependence of T -1 / 2 , while 
the RATE99 network assumes the rate coefficients to have 
no temperature dependence. More detailed studies show 



(Ncu feld et al .120051) . More work is clearly needed on these 
systems. 

At 10 K, 60% of the reactions present in both 
databases show a difference in rate coefficient lower than 
the considered uncertainty a condition which can be ex- 
pressed as 

hi/h2 



V2F, 



< 1, 



(1) 



where fci and k 2 are the rate coefficients from the different 
networks with k\ > k 2 and F is the uncertainty factor. In 
other words, almost half of the common reactions show a 
difference in rate coefficient not covered by the uncertain- 
ties in the rate coefficients, with 2% of them having an 
"extra" difference greater than three orders of magnitude 



V2Fi 



> 1000). These differences may not have strong 



consequences on the modeling results, however, if the re- 
actions concerned are not important or if formation and 
depletion reactions in one model are both larger to a sim- 
ilar extent than in the other. So, one must ask whether 
or not the differences between the abundances computed 
with the two databases are larger than the errors due to 
rate coefficient uncertainties. 

To answer this question, we ran Model 1 with the 
RATE99 network so as to compare with the osu.2003 re- 
sults. We obtained distributions of results for abundances 
that do not always overlap at the 2a level. In fact, eighty- 
four percent of the 373 common species show a disagree- 
ment at some times between the abundances computed 
with the two databases. No neutral species and only 62 
ions, such as S + , C^, CN + and Cg , overlap at all times. 
Note that the differences are due to difference in the val- 
ues of the rate coefficients and not in the differing reaction 
lists for the two networks, since we repeated the compari- 
son with a list of reactions common in both databases and 
obtained similar results. Fig. shows three examples: the 
CO, H 2 and CS molecules. The error in the abundance 
of CO at later times is very low because this molecule 
is the main reservoir of atomic carbon in the gas phase 
and both sets of reactions give the same results. At early 
times, between 6 x 10 3 and 2 x 10 5 yr, the abundance 
computed with osu.2003 is higher than the one computed 
with RATE99, reaching a ratio of 4 at 3 x 10 4 yr. The en- 
velopes defined by the error in the rate coefficients do not 
overlap. For H2O, the early and late abundances diverge 
by one order of magnitude and the envelopes do not quite 
overlap, whereas CS is produced to a much greater extent 
for the RATE99 case: the difference of up to 2 orders of 
magnitude at steady state far exceeds the envelopes of er- 
ror. Indeed, with RATE99, CS is found to be the reservoir 
of S: the difference between the 2 networks is thus quite 
profound and not only quantitative. 

We attempted to identify the reason for the discrep- 
an cy involving CS using the sensitivity method described 
in IWakelam et alJ l)2005|) . At 10 5 yr, of the 20 most im- 
portant reactions in the calculation of CS with the osu 
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rate coefficients by their values in RATE99, the CS result 
approaches the RATE99 abundance to within a factor of 
ps 3. Thus, one would argue that the sensitivity approach 
works well in deducing both the "important" reactions in 
the osu calculation and the reason for the difference with 
the results from the RATE99 calculations. But, one must 
be cautious because different results are obtained if one 
starts with RATE99: here, because of their smaller rate 
coefficients, ion-polar neutral reactions play a less signifi- 
cant role in the production and destruction of CS. Indeed, 
the reactions of importance for CS are quite different in 
the two networks. The introduction of the osu rates for 
a small number of " important" reactions in the RATE99 
network produces only a small change although the differ- 
ence in the rate coefficients of these reactions is not cov- 
ered by the uncertainty factor. The RATE99 abundance 
of CS is thus more stable against rate coefficient modifica- 
tions, which shows that a large number of reactions would 
have to be modified in order to change the result. 

For CO, which shows a factor of 3 difference at early 
times, we were not able to identify any main reactions 
responsible for this discrepancy. Once again, these small 
differences are then due to a large number of reactions 
with small variations of rate coefficients between the two 
databases. 

3.4. Comparison of uncertainty calculations 

IVasvunin et alJ l|2004f) did a similar study on dark cloud 
chemistry to ours using a previous version of the RATE99 
database with a flat distribution for the rate errors (see 
I Wakelam et al.l |2005, for a discussion concerning this 
point) and focused their results on the steady-state abun- 
dances. Although they defined the error in the abundance 
in a different way, we can compare their results with 
ours at 10 s yr. The authors used the full width at half 
maximum (FWHM) of the histograms of logX, which is 
~ 2^ A log X, as defined in Sect. [3 since most of the his- 
tograms have a Gaussian shap e (see the discussion about 
the definition of the error in IWakelam et ai1l2005F) . For 
mos t of the species , we f ound a lower uncertainty than 
did lVasvunin et all <j2004). As an example, the C4 mole- 
cule has a A log X of 1.45 at 10 8 yr in our study whereas 
IVasvunin et alJ found an FWHM of 2-3 for carbon clus- 
ters. The reason for this difference appears to be mainly 
due to the database used, because we found that the errors 
are generally higher with RATE99 than with osu. 2003, as 
can be seen in Fig. for the case of H 2 0. For the specific 
case of the molecule C4, we calculate that AlogX is 2 at 
10 s yr with RATE99, a value well above what we calcu- 
late with osu. 2003. The source of this difference in errors 
is unclear. 

4. Comparison with observations 

We have compared our theoretical results with some ob- 



Table 4. Observed abundances towards L134N. 



Species 


N(i)/N(R 2 ) 


Ref 


Species 


N(i)/N(R 2 ) 


Ref 


CH 


l(-8) 


(1) 


CN 


8.2(-10) 


(1) 


CO 


8.(-5) 


(1) 


CS 


1.7(-9) 


(2) 


NO 


6.(-8) 


(1) 


OH 


7.5(-8) 


(1) 


SO 


3.1(-9) 


(2) 


C 2 H 


<5.(-8) 


(1) 


C 2 S 


6.(-10) 


(1) 


H 2 S 


8.(-10) 


(1) 


HCN 


1.2(-8) 


(2) 


HNC 


4.7(-8) 


(2) 


OCS 


2. (-9) 


(1) 


SO2 


< 1.6(-9) 


(2) 


C 3 H 


3.(-10) 


(1) 


C 3 N 


< 2. (-10) 


(1) 


c 3 


<5.(-ll) 


(1) 


C 3 S 


< 2. (-10) 


(1) 


H2CO 


2. (-8) 


(1) 


H2CS 


6.(-10) 


(1) 


NH 3 


9.1(-8) 


(1) 


CH 2 CN 


< l-(-9) 


(1) 


C2H2O 


< 7.(-10) 


(1) 


C3H2 


2.(-9) 


(1) 


C 4 H 


l.(-9) 


(1) 


HCOOH 


3.e-10 


(1) 


HC 3 N 


8.7(-10) 


(2) 


CH 3 CN 


< l-(-9) 


(1) 


CH3OH 


3.7(-9) 


(2) 


CH3CHO 


6.(-10) 


(1) 


C2H3CN 


< l.(-io) 


(1) 


C3H4 


< 1.2(-9) 


(1) 


HC5N 


l.(-10) 


(1) 


HC 7 N 


2.(-ll) 


(1) 


HCO+ 


l.(-8) 


(2) 


HCS+ 


6.(-ll) 


(1) 


N 2 H+ 


6.8(-10) 


(2) 


H 2 CN+ 


< 3.1(-9) 


(1) 


H 2 


<3.(-7) 


(3) 


2 


< l-7(-7) 


(4) 


C° 


> l-(-6) 


(5) 









a a (-b) refers to a x 1 0"'' 
(11 Ohishi etld1(ll992h 
(2) [ Dickens et al.l ib oOO'l 
f3l| Snell et al.l JpOOCt) 
(^[ Paeani et alJ J2003j) 
(5) IStark et alJ fll99ftl 




log [t(yr)] log [t(yr)] 

Fig. 7. Theoretical abundance of OH (left) and CH 3 OH 
(right) as a function of time for Model 2. The dashed lines 
represent the observed values in L134N with an error of 
a factor of 3. The agreement between the observed and 
modeled abundance of OH is shown as a dark area whereas 
the disagreement for CH 3 OH is quantified by a distance 
of disagreement d (see Sect. 15.21) . 



are summarized i n Ta ble |H whil e thos e observed in TMC- 
1CP are listed in ISmith et ahl l)2004 Table 3). The goal 
of this comparison is to take into account both the uncer- 
tainty in the observed values and in the chemical model- 
ing. For the observed abundances, since all the abundances 
are not given with their corresponding errors, we assumed 
a standard uncertainty of ± a factor of 3 for the following 
reasons. First, telescope and atmospheric calibrations are 
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Table 5. Fraction of species reproduced by the different 
models and databases with the mean time of best agree- 
ment (expressed logarithmically). 



Factor 3 






Model 1 


Model 2 


Model 3 


L134N 


osu.2003 


73%, 5.7 


80%, 4.8 


78%, 4.7 




RATE99 


78%, 5.8 


83%, 5.9 


66%, 4.7 


TMC-1 


osu.2003 


57%, 4.4 


61%, 4.4 


76%, 4.9 




RATE99 


68%, 5.2 


70%, 5.7 


70%, 5.5 


Factor 5 






Model 1 


Model 2 


Model 3 


L134N 


osu.2003 


82%, 4.6 


87%, 4.7 


85%, 4.8 




RATE99 


83%, 4.8 


83%, 6.0 


74%, 5.2 


TMC-1 


osu.2003 


61%, 4.4 


67%, 4.4 


86%, 4.9 




RATE99 


72%, 5.6 


78%, 5.6 


76%, 5.3 



on the H2 column density, which is an additional source of 
uncertainty. The H2 column density is usually indirectly 
determined from the emission of other molecules via LTE, 
LVG, or Monte-Carlo models and usually different mole- 
cules give different estimates of H2 densities. Indeed, the 
emission of the molecules may come from different the 
layers of gas. In TMC-1CP for instance, it appears that 
C 18 may not come from the same volume of gas as the 
other molecules a nd the abundances compared to H2 may 
be overestimated (Pra tap et alJll997|) . Another reason for 
the high uncertainty of the observed values is that the 
inventory of the observed abundances come from several 
studies, in which several telescopes and approximations 
were used to compute the observed abundances, resulting 
in disagreements by at least a factor of three for some 
species. As a specific example, the S O abundance towards 
L134N (N position) is 3.1 x 10 ~ 9 inlDickens et alJ l|2000|) 
whereas it is 6 times higher in lOhishi et alj Tl992). Even 
though we think that an uncertainty of ± a factor of 3 
is reasonable, and consider it to be our standard observa- 
tional uncertainty, we also use a larger uncertainty of ± a 
factor of 5 to consider its effect. 

For the chemical model, we considered both the un- 
certainties in the rate coefficients and in the physical con- 
ditions, using Model 2 for both clouds. To compare the 
observed and modeled values, we first define agreement 
between the model and the observations to occur when the 
error bars of the calculated and observed values overlap. 
A more detailed approach is discussed below in Section 
15.21 An example of overlapping error bars at certain times 
only can be seen in Fig. 0for the case of OH in L134N, 
where the overlap is shown as a darker region. 

4.1. Results for L134N 

Figure |S] shows the times of agreement for each observed 
molecule towards L134N. With our standard uncertainty 
of a factor of three in the observed abundances, Model 2 



(ICS 
H 2 CO 



ocs 

HNC 
HCN 




log [t(yr)] 



log [t(yr)] 



Fig. 8. Agreement between the observed abundances to- 
wards L134N "N" peak) and the predictions of the chem- 
ical model (Model 2) as a function of time. Symbols in- 
dicate an agreement between the observational constraint 
and the model results. Stars are for the observed abun- 
dances whereas dots and diamonds are for the observa- 
tional upper and lower limits respectively. 



of agreement tapers off gradually at both younger and 
older ages. The peak number increases to 87% if an error 
of a factor of 5 is taken for the observed values. Model 1 
reproduces fewer species (73%), and the best agreement 
occurs later, at 5 x 10 5 yr (see Table |SJ). With Model 2, 
the molecules H 2 S, OCS, CH 3 OH, CH 3 CHO are under- 
estimated at all times. On the other hand, the observed 
abundances of HNC, C3H, C3O and NH 3 are reproduced 
by the model but not at the best range of ages. The three 
species C3H, HNC, and NH3 are in agreement with obser- 
vation for times very close to the best range, while C3O 
has not been detected and only an upper limit is known. 
The special cases of O2 and H2O, for which upper limits 
have been well studied, are discussed in Sect. 15.11 

4.2. Results for TMC-1 

Using Model 2, we can only reproduce at best 50% of 
the observed 52 molecular abundances in the CP peak of 
TMC-1. This percentage is even lower than the 67% ob- 
tained bv lSmith et alJ 1 20041) . who compared the results of 
the osu.2003 database with the same observations towards 
TMC-1 CP. These authors considered the modeling to be 
successful if the differences between the observed and the- 



oretic al values are less than a fa ctor of 10. |gmi^ h e^ al. 



2004JL iTerzieva fc Herbstl l(l998h and iRoberts fc Herbst 
2002T showed that by increasing the C/O elemental ra- 
tio, they were able to better reproduce the observations. 
We then ran a third model (Model 3, see Table identi- 
cal with Model 2 except that the initial abundance of O 
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C 2 H 2 
CHXN 



H 2 CS 
H z CO 



C 3 
C-N 



ocs 

HNC 
HCN 
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C 3 H 2 N + 
H 3 CN+ 
N 2 H 
HCS + 
HCO+ 
HC.N 
HC 7 N 
CH 3 C 4 H 
CH3C3N 
Oft. 
HC 5 N 

C 3 H 4 
C 2 H 3 CN 
CH,CHO 



C S H 
C^HZ 
CH 3 CN 
HNC, 
HCN 
HC 2 NC 
HCOOH 
C4H 



log [t(yr)] 



3 [t(yr>] 



Fig. 9. Agreement between the observed abundances to- 
wards TMC-1 ("CP" peak) and the predictions of the 
chemical model (Model 3) as a function of time. 



Here, seventy-six percent of the molecules are reproduced 
considering an observed uncertainty of a factor of 3 at a 
time of 8 x 10 4 yr. This number is increased to 86% at 
the same age for an observed uncertainty of a factor of 
5. The molecules OH, OCS, CH 3 OH, CH 3 CHO, C 2 H 3 CN 
and C3H4 are never reproduced by the model whereas the 
CN, CS, C 5 H, HC 9 N, C 6 H and CH3C3N abundances are 
not reproduced in the optimum age range. The observed 
abundance of OH in TMC-1 is quite uncertain since it 
was detected using a significantly larger beam than for 
the other species, and c ontamination from other regions 
of the cloud is probable IjOhishi et alJ ll9921. 

Looking at both clouds, we see that the abundance 
of the molecules OCS, CH 3 OH and CH 3 CHO are un- 
derestimated at all times by the model compared with 
the observed values. One explanation could be these 
species are formed on the grains and released in the ga s 
phase by non-thermal processes l|Markwick et al J 12000). 
Indeed, OCS and CH3OH are known to be present on 
the grain mantles with an abun dance (comparable to 
that of H 2 ) of - 10~ 7 for OCS llPalumbo et all Il997h 
and ~ 10~ 6 for CH3OH llChiar et al.lll996h . Work in 
progress by ICarrod et alJ l)submittedl ) with a gas-grain 
model and various non-thermal desorption mechanisms 
supports this view. The fact that H 2 S and NH3 are not 
reproduced correctly in L134N may have the same origin 
since H-enriched molecules are believed to form efficiently 
on grai ns although solid H 2 S has ne ver been detected on 
grains ijvan Dishoeck fc Blakel[l'99Sfl . For H-poor species 
such as CS, CN and HNC, some of the rate coefficients 
of the critical reactions forming and depleting them may 
be more uncertain than we have assumed, especially if the 
reactions have not been studied and the estimated uncer- 
tainties are hiding non-random errors. 



we find that with our sensitivity technique it is most of- 
ten difficult to isolate small numbers of very important 
reactions for those molecules with abundances we cannot 
reproduce. Rather, the general picture is that for many 
species, there are large numbers of reactions that are of 
relatively equal importance. 

There are some reactions deemed important by our 
sensitivity technique that are critical for all the species: 
these are the ionization of H 2 and He b y cos mic 
rays, as already noti ced by IVasvunin et alJ |2004h and 
IWakelam et alJ l|200fih . and raised in an oral presenta- 
tion by A. Markwick-Kemper in the 16th UCL Astronomy 
Colloquium (Windsor, UK, 2002). Ionization reactions in- 
volving cosmic rays, both direct and indirect, are, how- 
ever, better thought of as variable parameters rather than 
reactions with uncertain rates since it not in general the 
physics of the process but the flux of cosmic rays that is 
in doubt. 

4.3. RATE99 calculations 

Before definitively ascribing the disagreements to specific 
causes, it is useful to see how the comparison with ob- 
servations is affected by the use of the RATE99 network. 
Table |3] lists the overall level of agreement for both net- 
works. On balance, the results with RATE99 show slightly 
better agreement with observation and show it at later 
times. Unlike the osu.2003 results, the RATE99 results are 
in agreement for the saturated species methanol and, in 
some cases, acetaldehyde. Although this agreement would 
appear to weaken our argument that the saturated species 
are produced on grain surfaces, it should be noted that the 
methanol prediction of RATE99 is almost certainly wrong 
because the gas-phase synthesis of methanol used: 



cm 



H 2 



CH3OH; 



hi/, 



CH3OH+ + e™ — ► CH3OH + H, 



(2) 
(3) 



is assumed to be far too efficient for at least three rea- 
sons: (i) the radiative association reaction to produce the 
precursor ion CHaOH^ has b een measured to oc cur much 
more slowly than predicted l)Lucas et all [2002). (ii) the 
dissociative recombinatio n of this ion leads to methanol 
in only 6% of collisions l|Geppert et al.ll2006|) . and (iii) 
the water abundance used in the radiative association re- 
action is far greater than its measured upper limit. New 
calculations using the RATE99 network now agree with 
the osu.20 03 result that this sp ecies cannot be produced 
in the gas ijGeppert et al.ll200ql . 

The case of acetaldehyde may or may not be similar. 
This species is also thought to be produced via a radiative 
association 

H 3 0+ + C 2 H 2 — ► CH 3 CHOH+ + hi/, (4) 
followed by a dissociative recombination I Herbsdll987l) : 
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Fig. 10. Comparison between the theoretical abundances 
with error bars of H2O and O2 (Model 2) and the upper 
limits of the abundances found in L134N (see Table 



There is no direct evidence for either reaction, although 
the association does occur at hig h densities in the labora- 
tory via collisional stabilization (lHerbstlll987h . 

We conclude that there is substantial evidence that 
methanol is only produced on grains, although the case 
against gas-phase production of acetaldehyde is not quite 
proven. It is likely though that the high abundance of 
acetaldehyde in the RATE99 calculation stems at least 
partially from an overly slow destruction by ions since this 
species has a dipole moment. The percentages of species 
in agreement using the RATE99 database listed in Tabled 
do not include the agreement with these two species. 

5. Discussion 

5.1. Special cases of water and molecular oxygen 

An interesting result of this work is that it provides 
another manner in which the low upper limits to the 
abundances of gas-phase H2O and O2 in cold cores, de- 
tected by the SWAS and Odin satelli tes, can be explained 
l|Snell et al.ll20nnHPagani et al.ll2003h . Essentially, the up- 
per limit of the observed abundance (multiplied by a factor 
of three given our standard chosen observational uncer- 
tainties) must not be less than the lower error bar of the 
theoretical abundance. 

Let us first consider the case of L134N. In Fig. 
we plot the calculated Model 2 abundances for water and 
oxygen with error bars as functions of time and superim- 
pose the measured upper limits multiplied by a factor of 
three. The results show that both species are successfully 
reproduced in the optimum age range defined by the max- 
imum agreement with all species, although the case of O2 
is much more marginal since it is clearly overproduced at 
times after 10 5 yr. This result shows that interaction with 
grain surfaces is not necessarily required if we consider 
young clouds (< 10 5 yr). The abundance of O2 was found, 
however, to be lower t han the value for L 134N (< 10~ 7 ) 
in a dozen dark clouds l|Pagani et alJl2003ll . which we ex- 
pect to have different ages. It may then be necessary to 
invoke depletion processes to explain the low abun dance 
of O2 in quiescent cores fsee lRoberts &: Herbsti2002T) . The 
case of TMC-1 is singular since we used a model (Model 




Ot , , . 1 . . , 1 , , , ] 

2 4 6 8 




logMyr]) 

Fig. 11. Upper panel: the fraction of molecules repro- 
duced by the model as a function of time for both clouds 
(Model 2 for L134N and Model 3 for TMC-1). Middle 
panel: the distance D of disagreement between the ob- 
served and modeled abundances (see text for details). 
Lower panel: the fraction of reproduced molecules multi- 
plied by min(D) /D (minimum of the distance of disagree- 
ment divided by the distance of disagreement). We con- 
sider a factor of 3 uncertainty in the observed abundances. 



5.2. Age of the clouds: a more refined estimate 

Our models of quiescent cores are clearly not the complete 
picture since dynamical forces ar e producing and des troy- 
ing cores while chemistry occurs l|Garrod et alJl2005|) . For 
TMC-1 and L134N, the creation and chemistry have been 
occurring to some extent simultaneously. So, our pseudo- 
time-dependent calculations represent a crude approxi- 
mation, with frankly an uncertain zero-time condition. 
Even within this crude approximation, we have used a less 
than perfect criterion to determine the "chemical" age for 
TMC-1CP and L134N by simply maximizing the number 
of molecules with calculated abundances and error bars in 
agreement with observed values. Essentially, we have made 
no allowance for the quality of the agreement and the ex- 
tent of the disagreement for individual species. The results 
of this simple method are re-plotted in the top panel of 
Fig.[TT]for L134N (Model 2) and TMC-1CP (Model 3) as 
the fraction F of molecules in agreement vs time. In this 
representation, we can see that the best age for L134N is 
somewhat more distinct that for TMC-1 CP, but that the 
optimum ages lie in what previously was known as "early 
time" , well before steady-state conditions set in. 
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is defined, for the species not reproduced by the model, 
as the sum over all species of the distance di (in units of 
log A, see Fig.0) between the observational value and the 
theoretical one. For example, if for a particular species, i, 
X ts > X mo dei then the contribution to D is given by the 
expression 

di — l0gX o bs;min ~ log A mo del;max, (6) 

where the maximum and minimum values refer to the er- 
ror bars. The middle panel of Fig. II II shows the total dis- 
tance D as a function of time for the two clouds. This plot 
indicates that the disagreement with the model reaches a 
minimum, labeled min(D), somewhat after 10 4 yr. Finally, 
the lower panel of Fig. ^] shows a plot vs time of the ra- 
tio between the minimum distance, min(D), and D, mul- 
tiplied by the fraction of reproduced molecules, F. This 
quantity, which includes the strength of individual dis- 
agreements for specific species, tends to sharpen the best 
age, especially for L134N. In particular, we obtain a sharp 
maximum for L134N around 6 x 10 4 yr and a less marked 
maximum for TMC-1CP around 10 5 yr. These numbers 
are not strikingly different from a variety of other esti- 
mates, but must be taken with extreme caution. The ages 
derived analogously with the RATE99 set of reactions are 
somewhat larger, ~ 10 6 yr. 

6. Conclusions 

In this paper, we have reported the use of our Monte 
Carlo approach to uncertainties in calculated abundances 
to study the gas-phase chemistry of dark clouds, in partic- 
ular the well-studied sources TMC-1CP and L134N. We 
have utilized models in which the uncertainties in both 
rate coefficients and physical conditions are included; the 
latter can be thought of as either the actual observational 
uncertainties or physical heterogeneities in the sources. 
With a specific criterion for agreement between observa- 
tion and theory in which the errors in both techniques 
must overlap, we find that most but not species detected 
in the sources are reasonably well accounted for at so- 
called "early time," a far from novel result although one 
determined more rigorously than in previous approaches. 

A major goal of the study has been to gain some in- 
sight as to which failing of the simple picture utilized is 
the more important: the lack of grain-surface chemistry 
or the lack of a dynamical model of the sources. For sat- 
urated (H-rich) species, with reasonably understood syn- 
theses on grain surfaces, our results support earlier indi- 
cations that surface chemistry is an important, probably a 
critical omission, but for unsaturated species, the picture 
remains less clear because it is difficult to determine from 
our sensitivity analysis whether the discrepancies can be 
accounted for by poorly determined rates of critical chem- 
ical reactions. This difficulty stems from the fact that the 
chemistry of unsaturated species often seems to involve 



The fact that discrepancies for the calculated abun- 
dances of individual species between two chemical net- 
works used, our osu.2003 network and the RATE99 net- 
work, can be larger than their calculated random uncer- 
tainties indicates that more attention must be paid in 
the laboratory to specific classes of reactions, such as ion- 
polar neutral systems, before more definitive conclusions 
can be reached. Further complicating the issue are other 
concerns such as the proper initial make-up of the gas, 
what elemental abundances are reasonable, and the range 
over which the cosmic ray ionization rate £ can be varied. 
Nevertheless, our study indicates strongly that gas-phase 
models of static sources cannot represent the complete 
picture. More complex treatments, involving cycles of sur- 
face adsorption, reaction, and desorption, possibly coupled 
with dynamical histories of the sources, would appear to 
be required. 
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